# Making Comparisons 
sim_plant_dt_old =
  read.fst(
    path = here(paste0('Data/electricity-generation/output/quant-plant-dt.fst')),
    as.data.table = TRUE
  )
sim_plant_dt


plant_id_eia_check = unique(sim_plant_dt_old$plant_id_eia)[1235]

tmp = 
  merge(
    sim_plant_dt[plant_id_eia == plant_id_eia_check],
    sim_plant_dt_old[plant_id_eia == plant_id_eia_check], 
    by = c('quantile'), 
    all = TRUE
  )

tmp[,.(
  quantile, 
  damages.x, 
  damages.y,
  md_wecc.x, 
  md_wecc.y
)]


plant_gen_dt[
  plant_id_eia == plant_id_eia_check  &
  fit_net_gen_raw == fit_net_generation_mwh
]

merge(
  plant_gen_dt, 
  plant_dt[,.(plant_id_eia, namepcap)], 
  by = 'plant_id_eia'
)[fit_net_generation_mwh < namepcap]


sim_plant_dt[,.(
  min = min(fit_net_generation_mwh), 
  max = max(fit_net_generation_mwh)), 
  keyby = plant_id_eia
][,mean(min == max)]
sim_plant_dt_old[,.(
  min = min(fit_net_generation_mwh), 
  max = max(fit_net_generation_mwh)), 
  keyby = plant_id_eia
][,mean(min == max)]


# What is the marginal production of each fuel type by region?  
md_region_fuel_cat_sim_dt = 
  sim_plant_dt[,
    lapply(.SD, sum),
    keyby = .(quantile, tot_eload, nerc_adj, fuel_category),
    .SDcols = str_subset(colnames(sim_plant_dt),'md|mp')
  ] |> melt(
    id.vars = c('quantile','tot_eload','nerc_adj','fuel_category')
  ) %>%
  .[,.(
    nerc_adj_prod = nerc_adj,
    nerc_adj_demand = toupper(str_remove(variable, 'm(d|p)_')),
    fuel_category,
    variable = str_extract(variable, 'm(d|p)'),
    ic = fcase(
        toupper(str_remove(variable, 'm(d|p)_')) == 'TRE', 'Texas',
        toupper(str_remove(variable, 'm(d|p)_')) %in% c("CAL",'WECC'), 'West',
        default = 'East'
      ) |> factor(levels =c('West','Texas','East')),
    quantile, 
    tot_eload,
    value
  )] |>
  dcast(
    formula = nerc_adj_prod + nerc_adj_demand + fuel_category + ic + quantile + tot_eload ~ variable, 
    value.var = 'value'
  ) |>
  merge(
    quant_load_dt[
      square == FALSE & deviation == 'baseline', 
      .(tot_eload_ic = sum(excess_load)),
      by =.(
        quantile, 
        ic = fcase(
            toupper(nerc_adj) == 'TRE', 'Texas',
            toupper(nerc_adj) %in% c("CAL",'WECC'), 'West',
            default = 'East'
          ) |> factor(levels =c('West','Texas','East'))
        )
    ],
    by = c('ic','quantile')
  )


ggplot(
    md_region_fuel_cat_sim_dt[
      nerc_adj_prod == nerc_adj_demand
     & quantile %between% c(0.00,1.00)
    ], 
    aes(
      x = tot_eload_ic/1e3, y = md, 
      color = nerc_adj_demand, linetype = nerc_adj_demand
    )
  ) + 
  geom_line(size = 1.25) +
    scale_color_brewer(name = 'Region',palette = 'Dark2') + 
    scale_linetype_manual(
      name = 'Region',
      values = c('solid','twodash','dashed','longdash','solid','dashed','solid')
    )+
    labs(
      x = 'Excess Load (GWh)',
      y = 'Marginal Damage ($/MWh)'
    ) +
    facet_grid(cols = vars(ic), rows = vars(fuel_category), scales = 'free')+ 
    theme(
      legend.position = 'bottom',
      legend.key.width = unit(3, 'cm')
    )

# So it's a big drop in marginal damages at coal plants in RFC between 400 and 415 GWh
sim_plant_dt[
  nerc_adj == 'RFC' & 
  fuel_category == 'coal'
  #&tot_eload %between% c(400000,415000)
  ,
  .(
    n_plants = .N,
    n_plants_positive = sum(md_rfc > 0),
    V1 = sum(md_rfc)/sum(mp_rfc)), 
  keyby = quantile
] |>
ggplot(aes(x = quantile, y = V1)) +
geom_point()

|>
ggplot(
  aes(x = quantile, y = md_rfc)
) + 
geom_point()


sim_plant_dt[
  nerc_adj == 'RFC' & 
  fuel_category == 'coal'  & 
  quantile %in% c(0.92,0.93),.(
    quantile, 
    plant_id_eia, 
    md_rfc, 
    mp_rfc
  )
] |>
dcast(
  formula =  plant_id_eia ~ quantile,
  value.var = c('mp_rfc','md_rfc')
) %>%
.[,.(
  plant_id_eia, 
  mp_rfc_0.93, 
  mp_rfc_0.92,
  md_rfc_0.93, 
  md_rfc_0.92,
  diff_mp = mp_rfc_0.93 - mp_rfc_0.92,
  diff_md = md_rfc_0.93 - md_rfc_0.92
)] %>%
.[order(diff_mp)] %>%
ggplot(aes(x = diff_mp, y = diff_md)) + 
geom_point()

ggplot(
  aes(
    #x = quantile, 
    x = md_rfc, 
    color = as.character(quantile)
  )
) + 
geom_density()


# 2840 is problem plant
ggplot(
  sim_plant_dt[plant_id_eia %in% c('2840','6004')], 
  aes(x = quantile, y = md_rfc, color = plant_id_eia)
) + 
geom_line()

ggplot(
  sim_plant_dt_old[plant_id_eia %in% c('2840','6004')], 
  aes(x = quantile, y = md_rfc, color = plant_id_eia)
) + 
geom_line()

ggplot(
  sim_plant_dt[plant_id_eia %in% c('2840','6004')], 
  aes(x = quantile, y = mp_rfc, color = plant_id_eia)
) + 
geom_line()
ggplot(
  sim_plant_dt_old[plant_id_eia %in% c('2840','6004')], 
  aes(x = quantile, y = mp_rfc, color = plant_id_eia)
) + 
geom_line()

ggplot(
  sim_plant_dt[plant_id_eia %in% c('2840','6004')], 
  aes(x = quantile, y = md_rfc/mp_rfc, color = plant_id_eia)
) + 
geom_line()
ggplot(
  sim_plant_dt_old[plant_id_eia %in% c('2840','6004')], 
  aes(x = quantile, y = md_rfc/mp_rfc, color = plant_id_eia)
) + 
geom_line()


# Now have to figure out why production at that plant is decreasing so much 
prod_dt[,.(
  plant_id_eia, 
  id,
  square,
  est_extremum,
  excess_load,
  fit_net_gen_raw = fcase(
    est_extremum <= 0 | is.na(est_extremum), coef*excess_load, 
    est_extremum > 0 & excess_load <= est_extremum, coef*excess_load,
    est_extremum > 0 & excess_load > est_extremum, coef*est_extremum,
    default = NA
  ), 
  coef*excess_load, 
  coef*est_extremum
)][!is.na(est_extremum)]

# So where does est extremum come from? 
model_dt_long[plant_id_eia == plant_id_eia_in]
0.00506232804274/(2*0.00000002203211)
